function [t,y] = euler(odefun,tspan,y0,h)
    t = tspan(1) : h : tspan(2);
    y(1) = y0;
    for i = 1:length(t) - 1
        y(i + 1) = y(i) + h *odefun(t(i),y(i));
    end

end